function errors

%  Plots the error as a function of the number of grid points 
%  for the BVP
%       y'' + p(x)y' + q(x)y= f(x)   for xL < x < xr  '
%  where 
%      y(xl) = yL  and y(xR) = yR

%  p=0, q=-1, f=sin(2*pi*x)
%  xL=0, yL=0  and  xR=1, yR=0

% clear all previous variables and plots
clear *
clf

% set boundary conditions
	xL=0; yL=0;
	xR=1; yR=0;

%  x-location where error is calculated
x3=1/3;  exact=-sin(2*pi*x3)/(1+4*pi*pi);
%x4=1/4; exact=-sin(2*pi*x4)/(1+4*pi*pi);

% iteration to determine the error as number of points is increased
ii=0;
for i=0:3

	% set number of points along the x-axis
	n=3*10^i+1;
	%n=4*10^i+1;
	ii=ii+1; 
	points(ii)=n;

	% index for x=x(ixpoint) where error is to be measured
	ixpoint=10^i+1;

	% generate the points along the x-axis, x(1)=xL and x(n)=xR
	x=linspace(xL,xR,n);
	x(ixpoint)
	h=x(2)-x(1);

	% calculate the coefficients of finite difference equation
	a=zeros(1,n-2); b=zeros(1,n-2); c=zeros(1,n-2); z=zeros(1,n-2);
	for i=1:n-2
		a(i)=-2+h*h*q(x(i+1));
		b(i)=1-0.5*h*p(x(i+1));
		c(i)=1+0.5*h*p(x(i+1));
		z(i)=h*h*f(x(i+1));
	end;
	z(1)=z(1)-yL*b(1);
	z(n-2)=z(n-2)-yR*c(n-2);

	% solve the tri-diagonal matrix problem
	y=tri(a,b,c,z);
	y=[yL, y, yR];

	error(ii)=abs(exact-y(ixpoint));
	ye=-sin(2*pi*x)/(1+4*pi*pi);
	errorM(ii)=norm(ye-y,inf);
end;

loglog(points,errorM,'--b*','LineWidth',1,'MarkerSize',7)
hold on
loglog(points,error,'--sr','LineWidth',1,'MarkerSize',7)
grid on
set(gca,'MinorGridLineStyle','none')

% define title and axes used in plot
xlabel('Number of Grid Points','FontSize',14,'FontWeight','bold')
ylabel('Error','FontSize',14,'FontWeight','bold')
legend(' Max Error',' Error at x = 1/3',1);
title('BVP: Error in Example 1','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  

hold off


function g=q(x)
g=-1;

function g=p(x)
g=0;

function g=f(x)
g=sin(2*pi*x);


% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end